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We review c-field methods for simulating the non-equilibrium dynamics of degenerate Bose gases beyond 
the mean-field Gross-Pitaevskii approximation. We describe three separate approaches that utilise similar nu- 
merical methods, but have distinct regimes of validity. Systems at finite temperature can be treated with either 
the closed-system projected Gross-Pitaevskii equation (PGPE), or the open-system stochastic projected Gross- 
Pitaevskii equation (SPGPE). These are both applicable in quantum degenerate regimes in which thermal fluc- 
tuations are significant. At low or zero temperature, the truncated Wigner projected Gross-Pitaevskii equation 
(TWPGPE) allows for the simulation of systems in which spontaneous collision processes seeded by quantum 
fluctuations are important. We describe the regimes of validity of each of these methods, and discuss their rela- 
tionships to one another, and to other simulation techniques for the dynamics of Bose gases. The utility of the 
SPGPE formalism in modelling non-equilibrium Bose gases is illustrated by its application to the dynamics of 
spontaneous vortex formation in the growth of a Bose-Einstein condensate. 



I. INTRODUCTION 

It may seem counter-intuitive that the dynamics of Bose- 
Einstein condensation, a phenomenon that occurs due to 
quantum statistics, can often be well-described by an equa- 
tion of motion for a classical field. However, there are in fact 
two distinct classical regimes [1| that can potentially be con- 
fused. The first is the classical-particle regime — where the 
de Broglie wavelength of the atoms in a gas is much smaller 
than the typical interparticle spacing, and the particles can be 
treated as 'billiard balls'. The second is the classical-field 
regime, in which the modes of a bosonic quantum field are 
highly occupied ((Nk) » 1), such that the discrete set of inte- 
ger mode occupations can be approximated by a continuum. 
This occurs near quantum degeneracy for the ultra-cold Bose 
gas, and in these circumstances a description using a classical 
wave equation may be appropriate. There are two cases worth 
distinguishing: The Gross-Pitaevskii equation (GPE) J2] |3) 
for a pure Bose-Einstein condensate (BEC) at zero tempera- 
ture approximates the highly-occupied condensate as a clas- 
sical field, analogously to the representation of the quantum 
field of laser light by a classical electromagnetic wave ||4). 
However, just as the long-wavelength limit of the Planck dis- 
tribution of blackbody radiation — in which the distinction 
between individual field quanta is unnecessary — can be de- 
scribed by the classical Rayleigh-Jeans law [5 1, the GPE can 
also provide, within appropriate limits, a description of the 
finite-temperature atomic Bose field. This Chapter describes 
methodologies in which an equation of motion for a classical 
field provides a beyond-mean-field description of the degen- 
erate Bose gas. 

A number of authors have proposed that a classical-field 
treatment of the Bose gas is appropriate to describe the dy- 
namics of condensate formation. Near degeneracy, a quantum 
Boltzmann equation description becomes inadequate due to 
the growing coherences between the modes. Kagan, Svis- 
tunov, and Shlyapnikov suggested that when the inequality 
Nk = {a^tk) » 1 is satisfied for a large number of the low- 
energy modes of the gas, the Bose field is more appropriately 
described by a Gross-Pitaevskii equation ll6l fT0l . Other au- 



thors, including Stoof [11 1 and Gardiner et al. |[T2l[T3l . came 
to realise the need for a classical-field description from differ- 
ent starting points. 

The first numerical simulations of the GPE demonstrating 
the kinetics of Bose condensation in a classical-field model 
were performed by Damle et al. Ifl4l . This was followed by 
work by Marshall et al. fl3], Stoof et al. [16] Q31, Goral et 
al. EH, Sinatra et al. (HI, Davis et al. ll20H22l . and Berloff 
and Svistunov ||23l , who all simulated the Bose gas at finite 
temperature using methods based on the GPE. Collectively, 
this type of approach has become known as the 'classical- 
field method' ifTTl l24l l25l . While the simplest variants of this 
method entirely neglect the effects of quantum fluctuations, 
they have the significant advantage that they treat the classi- 
cal fluctuations of the field non-perturbatively, and hence can 
be applied in the critical regime around the Bose-condensation 
phase transition 1 26-3 1 1 and the fluctuation regime of the two- 
dimensional Bose gas Il32l - l35l . The related truncated Wigner 
methods [24, 36-39 ) incorporate the leading-order effects of 
quantum fluctuations into a GPE-like description [40|. 

In this Chapter we briefly outline the derivation of the 
stochastic projected Gross-Pitaevskii equation (SPGPE) for- 
malism — for a recent, more extensive review see Ref. 11241 . 
The SPGPE is an equation of motion for a classical field 
(c-field) describing the low-energy, highly-occupied modes of 
a Bose gas, coupled to a bath of high-energy atoms that is as- 
sumed to be close to thermal equilibrium. It offers a powerful 
framework for the study of both equilibrium correlations — 
including condensation, anomalous correlations, and critical 
fluctuations — and non-equilibrium dynamics of the finite- 
temperature Bose gas. We also discuss two other related meth- 
ods: The projected Gross-Pitaevskii equation (PGPE) formal- 
ism is obtained from the SPGPE upon neglecting the coupling 
of the c-field to the bath, and the projected truncated- Wigner 
method (TWPGPE) incorporates the leading-order effects of 
quantum fluctuations l40l at low temperatures, where such 
fluctuations may be important. We describe each of these 
methods and their regimes of validity, and discuss their re- 
lationship to the work of other authors. Finally we describe 
the application of the SPGPE to the modelling of spontaneous 
vortex formation in condensate growth. 
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II. METHODOLOGY 



A. Outline of Derivation 



The classical Rayleigh-Jeans law for thermal black-body 
radiation is a good approximation for the long-wavelength 
modes of the electromagnetic field, but results in the well- 
known ultra-violet catastrophe at short wavelengths Q • This 
suggests that a classical treatment may provide a good approx- 
imate description of the atomic Bose field when restricted to 
the low-energy field modes. We therefore divide the Bose field 
operator *P(r) into a low-energy part that will be treated clas- 
sically (the coherent or c-field region C), and a high-energy 
part that will be treated quantum mechanically (the incoher- 
ent region I) 



T(r) = T c (r) + Y,(r). 



(1) 



This division is effected by defining a projection operator onto 
the c-field region 

W(r)} = J] fn(r) f dr' <p* n (r')f(r'), (2) 



such that 



*c(r) = !Pcmr)} = 2^" (r) ' 



(3) 



neC 



where the <p„(r) are the eigenvectors of the appropriate single- 
particle Hamiltonian, with eigenvalues s„. The division be- 
tween the C and I regions is made at a cutoff energy E cut , 
which is chosen such that the modes in C are classical, that is 
<jVjt> » 1 for the highest energy mode in C. The particular 
choice of E cut will be dependent on the thermodynamic pa- 
rameters of the initial state of the system, but we require that 
the physical observables we calculate should be insensitive to 
the exact value chosen for the cutoff 1 30 1 . 

Our goal is to derive a computationally tractable equation of 
motion for the c-field region. To achieve this we make use of 
the methods of open quantum systems BP . and treat the inco- 
herent region I as a thermal and diffusive reservoir (or bath) to 
which the field ^(r) is coupled [ 13 1. We assume that the bath 
density operator is thermal and quantum Gaussian [41 1, so that 
all operator products factorise into products of one-body cor- 
relation functions. Furthermore, we assume that these corre- 
lations are well represented by a semiclassical description in 
terms of a one-body Wigner function F\(r, K) RJ2). In princi- 
ple, an equation of motion for F\(r, K) could be derived using 
the methods of kinetic theory, and this is a topic of current 
research. Here we assume thermal equilibrium such that 



^i(r,K) 



\l/{expmT,K)-n)/k B T] 
10 



1} g(r,K)>£ cut , 
e(r, K) < Ecat, 
(4) 

where s{r, K) is the semiclassical energy 11431 of a particle 
with position r and wavevector K. The I region is then com- 
pletely characterised by its chemical potential // and temper- 



ature T (and any other appropriate thermodynamic Lagrange 
multipliers [44 J), and these parameters will arise in the equa- 
tion of motion for the classical field. 

A master equation for the coherent-region density op- 
erator can be obtained by tracing out the incoherent re- 
gion, using standard methods from quantum optics BP . A 
'high-temperature' approximation is then made for the mas- 
ter equation, requiring that all the eigenfrequencies of the 
C-region evolution are small compared to the temperature; 
i.e., H\ej\/knT <K 1 [13] (in practice we find h\ej\ < fi + 
E C ut l45l ). The master equation can then be mapped to a 
generalised Fokker-Planck equation for the Wigner quasi- 
probability distribution using the operator correspondences 
described in Refs. 11241 14TI . The binary interaction term in 
the standard Bose-field Hamiltonian leads to terms involving 
third-order derivatives of the Wigner distribution with respect 
to the phase-space variables (36), so there is no exact repre- 
sentation of the resulting Fokker-Planck equation in terms of 
a stochastic differential equation (SDE). However, for highly 
occupied c-field modes we may make the truncated Wigner 
approximation (TWA), in which the third-order derivatives 
are neglected l36ll . The resulting (Stratonovich form J46)) 
stochastic differential equation for the classical-field (D(r) is 
the SPGPE 



(5a) 
(5b) 



dO(r)=P c \-^H O(r)dt 



G(r) 
k*T 



(jx - Ho)$>(r) dt + dW G (r, t) 



+ I dr'M(r - r') ^—<&(r)dt + id>(r)dW M (r, t) 



k B T 



where 



2m 



V 2 + Vext(r)+^(r)| 2 



(5c) 



(6) 



is the Gross-Pitaevskii operator, with g the three-dimensional 
interaction strength. All terms on the RHS of the SPGPE in- 
volve the c-field projection operator Pq, which formally re- 
stricts the dynamics of cD(r) to the region C. Aside from the 



projection, the first line of the SPGPE, Eq. ( 5a i, is the standard 
zero-temperature GPE. 

Equations (|5j3)-(j5j;) describe the coupling of the c-field to 
the I region. Equation pb| corresponds to growth processes, 
in which two I-region atoms collide and one is scattered into 
C, and the corresponding time-reversed processes. The local 
c-field amplitude grows if the value of Ho<l>(r) is smaller than 
/LiO(r), and vice versa. The growth rate G(r) is given by a col- 
lision integral over the I region — see Eq. (171) of Ref. |2"4"1 . 



Equation ( 5c I corresponds to scattering processes in which 
a particle in the I region scatters off one in the C region, leav- 
ing the total number in each region unchanged. These pro- 
cesses couple the bath to the divergence of the c-field current 



jc(r) = ^([V0>*(r)]O(r) - <D*(r)V(D(r)), 



2m 



(7) 
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and are characterised by the scattering amplitude M(r - r'), 
which is given by another collision integral [Eq. (174) of 
Ref. El]. 

The SPGPE involves two explicit stochastic terms: a com- 
plex noise process dWc{r, t) associated with growth, and a 
real noise process cIWm(t, associated with scattering. The 
non-zero correlations of these Wiener increments are 

(dW* c (r, t)dW G (r', f)> =2G(r)(5 c (r, r')dt, (8) 
(dW M (r, t)dW M (r', t)) =2M(r - r')dt, (9) 

where 5c{r, r') = 2„ec tfin( r )'P*r,( r ') is the kernel of the projec- 
tion operator Vq and acts as a Dirac delta function within the 
c-field region. 

As the SPGPE is derived within the TWA, the initial con- 
dition 0(r, t — 0) should in principle include a representa- 
tion of quantum fluctuations; i.e., half a particle of noise per 
mode |24|. Distinct samples of this initial 'vacuum noise' and 
the dynamical noise processes generate distinct SPGPE tra- 
jectories, and in the Wigner formalism moments of the c-field 
O(r) calculated from an ensemble of such trajectories corre- 
spond to quantum-statistical expectation values of symmet- 
rically ordered field-operator products [36|. In practice, the 
evolution of the field is rapidly dominated by classical fluctua- 
tions which overwhelm the contributions of the vacuum noise, 
and the formal operator-ordering corrections of the Wigner in- 
terpretation are therefore usually neglected; i.e., moments of 
the c-field are interpreted directly as moments of the quantum 
field, which is of course a reasonable approximation in the 
classical-field limit in which the SPGPE is derived. Moreover, 
individual trajectories are often interpreted as 'typical' reali- 
sations of a particular experiment — see Ref. [24] for further 
discussion. 

Although fundamentally describing an open dynamical sys- 
tem, the SPGPE can be used to non-perturbatively calculate 
the equilibrium properties of finite-temperature Bose gases. 
Once any transient dynamics have subsided, the SPGPE tra- 
jectories sample the grand-canonical equilibrium of the sys- 
tem in the classical-field limit, in the spirit of the well-known 
Langevin equation approach to sampling thermal distribu- 
tions flTl . and have been utilised, for example, to calculate 
the thermodynamic properties of trapped ID Bose gases 11481 . 
Dynamical calculations using the SPGPE method have mostly 
focussed on the decay of vortices |49 50] and the formation 
of condensates 144, 5T1. 



B. Numerical Implementation of the SPGPE 

Applications of the SPGPE have thus far all focussed on a 
numerically expedient simplification of the full formalism — 
the 'simple growth SPGPE' — which utilises two further ap- 
proximations: First, the growth rate G(x) is basically constant 
over the bulk of the c-field region [44 1, so we neglect its spatial 
variation. Second, the scattering terms appearing in Eq. ( 5c i 
are neglected, on the basis that these terms are small unless 
large c-field currents are present, which does not occur near 
equilibrium. 



Numerically, the non-local nature of the scattering term in- 
troduces significant complexity, and as the associated noise 
is multiplicative, first- or second-order (and hence inefficient) 
integration algorithms would be required for stochastic con- 
vergence. The simple growth SPGPE is, by contrast, an SDE 
with additive noise, which are typically more stable than SDEs 
with multiplicative noise [52]. We have integrated the simple 
growth SPGPE using an interaction-picture method ll46ll . and 
a fourth/fifth-order Runge-Kutta algorithm with an adaptive 
step size ll53l . Such an algorithm is not strictly convergent 
for SDEs, but in the limit of weak noise it is an acceptable 
approximation 1 52 1 . 

The implementation of the SPGPE is greatly simplified if 
the c-field <D(r) is expanded over the single-particle basis in 
which the projection operator Vc is defined. In order for 
the cutoff defined in terms of the basis {yj„(r)} to be sensi- 
ble, the highest-energy excitations in the C region should be 
single-particle-like (i.e., relatively unaffected by interactions). 
Fortunately, for many situations of interest this condition can 
be satisfied simultaneously with the classical-field condition 
(N k ) » 1 El 1231 EH- 

Many of the terms in the SPGPE are spatially local, and 
so efficient transformations between the modal representation 
and position space are required. For homogeneous systems a 
plane-wave basis is appropriate, and the necessary transforma- 
tions can be implemented using fast Fourier transforms ETI . 
For harmonically trapped systems it is possible to use Gaus- 
sian quadrature methods to effect numerically exact transfor- 
mations Il22l l55l . Although the computational cost of these 
transformations scales (in 3D) as M 4 (where M is the num- 
ber of modes in each dimension) as opposed to the M 3 log M 
scaling of the fast Fourier transform 11551 [56], the value of 
M required is somewhat smaller than is needed to accurately 
represent the same harmonic oscillator states on a Cartesian 
grid El. 



C. Projected Gross-Pitaevskii Equation (PGPE) 

The projected Gross-Pitaevskii equation (PGPE) is a sim- 
plification of the SPGPE obtained by retaining the high- 
occupation validity condition (Nk) » 1 for the modes of the c- 
field region, but neglecting their coupling to the I region. This 
results in a closed system, and the resulting Hamiltonian evo- 
lution of the field <I>(r) conserves the energy, normalisation, 
and any other first integrals which may be present, such as the 
momentum l20l I2D . angular momentum 11581 . or spinor-gas 
magnetisation (59]. 

The PGPE has mostly been used for calculating the equi- 
librium properties of Bose gases at finite temperature 11201 - 
|22l H9H2D El EH l58ti6ll. In practice, the PGPE seems to 
generate a suitable ergodic time evolution; i.e., given an ar- 
bitrary initial state, the trajectory of the field densely covers 
the microcanonical ensemble defined by the conserved en- 
ergy and other first integrals l62l . Once equilibrium is estab- 
lished, time averages of the field can be substituted for ther- 
modynamic ensemble averages lfT8ll20ll2T1l . Rugh has shown 
generically [63] that derivatives of the classical microcanoni- 
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cal entropy can be reformulated as averages of functions of the 
phase-space variables over the appropriate constant-energy 
surface. This methodology has been applied to the PGPE, 
and allows for the rigorous calculation of the temperature and 
chemical potential of the system by time averaging ll2"9ll64"l . 
The above-cutoff region can then be described in a mean-field 
approximation [ 30 1 . 

Because the PGPE neglects collisional processes that trans- 
fer population between the coherent and incoherent regions, 
it is likely to underestimate damping rates (for example), and 
the PGPE dynamics are potentially sensitive to the value of the 
cutoff [65|. However, dynamical calculations within a pure 
PGPE formalism are able to provide useful insights into the 
dynamics of degenerate Bose-gas systems in situations where 
a precise identification of the method with the full field theory 
is impractical [45 58 66 1. The PGPE has also been used to es- 
tablish the connection between c-field methods and more tra- 
ditional theoretical methods based on U(l) symmetry break- 
ing ED mi. 



D. Truncated Wigner Projected Gross-Pi taevskii Equation 
(TWPGPE) 

The truncated Wigner projected Gross-Pitaevskii equation 
(TWPGPE) 1 24 1 is an implementation of the truncated Wigner 
approximation [19, 36 , 37 1 that is computationally identical to 
the PGPE, but has a distinct regime of validity. Whereas the c- 
field in the SPGPE and PGPE methods constitutes an approxi- 
mation to the quantum Bose field in the sense that it describes 
only highly occupied field modes, the TWPGPE includes a 
representation of the quantum fluctuations, and can thus de- 
scribe weakly occupied modes as well. However, the formal 
Wigner interpretation of the trajectories quickly becomes in- 
validated at the high temperatures at which the SPGPE and 
PGPE are applicable. 

The TWPGPE is distinguished from the truncated Wigner 
method described in Refs. 1 19, 36, 37] by the implementation 
of a high-energy cutoff using a projection operator. This pro- 
jection operator formalises the coarse-grained effective-field- 
theory description of the Bose field |28 | underlying the de- 
scription of interatomic interactions by a 'contact' potential 
(see, e.g., Ref. [68]), and provides a rigorous basis for the 
addition of a finite density of quantum fluctuations to what 
would otherwise be a divergent local field theory |69l . An ap- 
propriate choice of projection operator can mitigate some of 
the spurious effects of quantum noise that may arise within the 
TWA in two and three dimensions, by effectively reducing the 
number of simulated modes (and hence the amount of noise) 
significantly 11701 . 



III. VALIDITY ISSUES 

A well-established validity condition for the TWA at T = 
is that the number of particles being simulated should be 
somewhat larger than the number of modes in the c-field (37]. 
Intuitively, the real population of each mode should dominate 



the half-quantum vacuum occupation of the mode, in order 
for the neglect of the quantum processes represented by third- 
order derivatives to be valid. This has been formalised by 
Polkovnikov BUI , who developed a perturbation expansion 
for quantum dynamics around the classical-field limit, and 
found that the TWA is obtained as the first-order correction 
to the classical (GPE) dynamics. Norrie et al. [69 1 argued 
that fundamentally, the local density of real particles should 
be large compared to the density of vacuum fluctuations. As 
the terms neglected in the TWA are only significant where 
the total particle density is large, the inaccuracy of the evo- 
lution in regions of low real-particle density is therefore of 
comparatively little consequence for the overall accuracy of 
the method. 

Over time the error associated with the Wigner truncation 
grows, and at long times, the ergodic character of the PGPE 
causes its solutions to thermalise to a classical microcanoni- 
cal equilibrium, so that a formal Wigner interpretation of the 
trajectories is no longer available. Sinatra et al. noted 11371 
that thermalisation of the initial noise population will lead 
to a spuriously high equilibrium temperature, and the rates 
of damping processes in the field may therefore be overesti- 
mated. However, in non-equilibrium scenarios, this effect may 
be of little consequence, even in simulations starting from zero 
temperature B31 l57l . In the limit of classically thermalised 
fields such as those described by PGPE and SPGPE equilib- 
ria, in which the formal Wigner moment correspondences are 
neglected, the fundamental validity condition is the mode oc- 
cupation condition (Nk) » 1. In practice, this is well satisfied 
at temperatures ranging from just above the critical tempera- 
ture T c down to about 0.5 T c 1711 . 



A. Relevance to Other Theories 

The SPGPE method is closely related to the stochastic GPE 
derived by Stoof [ 1 1 1 within a Keldysh path-integral formal- 
ism. In Stoof's approach the effects of dissipation, which are 
described using quantum-optical tools in the derivation of the 
SPGPE, appear as self-energy terms (cf. Ref. [72 1). Proukakis 
and Jackson [73] have shown that the Stoof SGPE is func- 
tionally equivalent to the SPGPE, when the scattering term 
and projection operator are neglected. However, the precise 
relationship between the two formulations remains somewhat 
unclear, as it does not appear that a strict division between 
the classical modes described by the Stoof SGPE and the 
modes that form the 'reservoir' to which the classical field are 
coupled is made in the Stoof approach. Despite these differ- 
ences, comparisons of numerical results of the simple growth 
SPGPE and the Stoof SGPE show close agreement in equilib- 
rium G3. 

Ergodic classical-field simulations for equilibrium correla- 
tions, similar to the PGPE, have been considered by other au- 
thors lf]~8l |25~1 rzST - FTTl . These simulations typically do not in- 
clude an explicit projector, leaving the cutoff to be determined 
by the spatial resolution of the numerical grid. However, cal- 
culating the nonlinear term of the GPE on such a grid may lead 
to numerical aliasing — a well-known issue in computational 
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FIG. 1: Spontaneous vortices in the formation of a Bose-Einstein condensate, a. Squares: number of condensate atoms as a function of time 
measured in the experiment. Solid line: Condensate number from SPGPE simulations with parameters: T\ = 45 nK, T f = 34 nK, ^ = %u r , 
Hi = 25na) r , TiG(r)lk^T = 0.005. Dashed line: probability of finding one or more vortices in the condensate as a function of time, averaged 
over 298 trajectories. The shaded area indicates the statistical uncertainty in the experimentally measured vortex probability at t = 6.0 s. It 
was observed in experiment that there was no discernible vortex decay between 3.5 s and 6.0 s b. Experimental absorption images taken after 
59 ms time-of-flight showing the presence of vortices, c. Simulated in-trap column densities at t = 3.5 s (indicated by the vertical dotted line 
in a.) d. Phase images through the z = plane, with plusses (open circles) representing vortices with positive (negative) circulation. Adapted 
with permission from C. N. Weiler et al, Spontaneous vortices in the formation of Bose-Einstein condensates, Nature 455, 948 (2008) |51 1. 



physics ll78l — particularly as the classical-field approxima- 
tion requires that all simulated modes are significantly occu- 
pied. These authors also employ somewhat ad hoc refine- 
ments to estimate Bose-field correlations from a classical- 
field model. For example, Sinatra and co-workers f76l 1771 
have suggested modifying the single-particle Hamiltonian in 
the interacting classical-field equation of motion, such that the 
equipartition of system energy would produce a bosonic num- 
ber distribution in the ideal-gas limit. Similarly, Brewczyk, 
Gajda, Rzazewski and co-workers [25] have attempted to de- 
fine an 'optimal momentum cutoff', by adjusting the com- 
putational grid spacing in their simulations of the interacting 
Bose gas, such that the occupations of the condensate or other 
modes match those of the ideal gas at the same temperature 
11791 [80l . In our view, the division of the Bose field into clas- 
sical and incoherent regions, implemented by means of a pro- 
jection operator that properly restricts the classical-field ap- 
proximation to modes of the field for which it is valid, and 
allows for a quantitative mean-field description of the com- 
plementary high-energy part of the field, is the only classical- 
field approach which facilitates truly quantitative descriptions 
of experimental systems (see, e.g., Ref. |30|). 

The ZNG method for finite-temperature Bose-Einstein 
condensates described by Zaremba, Nikuni, and Griffin llSTI 
bears some superficial similarity to the SPGPE approach, in 
that it involves a generalised Gross-Pitaevskii equation which 
features damping terms representing the effects of the thermal 
cloud. However, the ZNG approach is based on the funda- 
mentally different perspective of assuming the existence of 
a well-defined condensate mode, which is introduced via a 
symmetry -breaking ansatz 1821 . By contrast, in the SPGPE 
approach no such assumption of condensation is made, and 
the condensate (if one is present) must be extracted a pos- 
teriori from the field correlations 11221 [83l . Further discus- 
sion of the different philosophies of the c-field and ZNG ap- 



proaches can be found in Refs. 1 84- 1 and [85]. The main advan- 
tage that the ZNG method offers is a description of the cou- 
pled dynamics of the condensate and the full thermal cloud, 
in which the kinetics of the latter are modelled in a Boltz- 
mann equation approach. Such a coupled condensate-cloud 
description appears to be essential for accurately describing 
certain collective oscillations of the gas at high temperatures 
ll65l l86l 1571 . However, as it is based on the assumption that a 
well-defined condensate exists, the ZNG approach is not ap- 
plicable to more general scenarios involving low-dimensional 
systems Il34ll48l[88l[89l , regimes of turbulent matter-wave dy- 
namics Ir23ll45l l90l, or non-equilibrium passage through the 
transition to condensation lfT6ll44ll5Tl . 



IV. APPLICATION 

The first experiments on 3D condensate formation observed 
BECs to grow essentially quasistatically as smooth Gross- 
Pitaevskii-like states I9T1 1921 . and were successfully mod- 
elled using kinetic-theory approaches [93-95 1. Later exper- 
iments by Shvarchuck et al. l96l and Hugbart et al. 11971 stud- 
ied the formation dynamics of non-equilibrium quasiconden- 
sates. More recently, the Anderson group at the University 
of Arizona evaporatively cooled Bose gases rapidly from near 
degeneracy to BEC in a weakly oblate harmonic trap llBTl . 
When imaging along the vertical (z) direction, they found that 
approximately 25% of BEC column-density images featured 
holes consistent with vortices in the condensate. 

Simulating the Anderson group's experiments with the sim- 
ple growth SPGPE has helped elucidate the phenomenon of 
spontaneous vortex formation [51 1. The initial states for the 
simulations were sampled from an equilibrium ensemble near 
degeneracy, with a temperature T and chemical potential /j. 
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Evaporative cooling was modelled by suddenly changing the 
bath parameters T and yt to experimental values measured im- 
mediately after the most rapid phase of condensate growth. 
The dimensionless growth rate TiG{r)lk^T was set to a con- 
stant in time and space chosen so that the simulated conden- 
sate growth curve matched that observed experimentally, as 
shown in Fig. [T] A total of 298 trajectories were simulated, 
with each interpreted as a separate experimental run. 

The simulated probability of observing a single vortex in 
the resulting condensate agrees well with the experimental re- 
sults, as indicated in Fig. [T] The majority of spontaneously 
formed vortices aligned with the z axis, and hence were eas- 
ily observed in the column-density images. The vortices were 
found to live for many seconds in the experiment, somewhat 
longer than in the simulations. However, the temperatures ul- 
timately reached in the experiment are below those for which 
the SPGPE is applicable, and so this additional cooling af- 
ter the initial growth phase was not simulated. The discrep- 
ancy between experimentally measured vortex lifetimes and 
the SPGPE simulation results is therefore consistent with the 
dependence of vortex decay rates on the system tempera- 
ture l|49ll98l[T00l . These results demonstrate that the SPGPE 
is a powerful tool for the quantitative modelling of the non- 
equilibrium dynamics of Bose gases at finite temperature. 



V. RELEVANCE TO OTHER SYSTEMS 

The SPGPE could potentially be useful for simulating non- 
equilibrium exciton-polariton systems, which are only ever in 
steady state due to the pumping and finite lifetime of the quasi- 
particles (see, e.g., Ref. [ 101']). However, a more general rep- 
resentation of the I region would likely be necessary, as the 
reservoir will not in general be in equilibrium, and the scat- 
tering term neglected in the simple growth SPGPE could play 
an important role in establishing the steady state of the C re- 
gion. Wouters and Savona have applied methodology similar 
to the SPGPE 1 102 1 to investigate superfluidity in polariton 
systems 11 1 031 - 
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